﻿#pragma once
#include "BaseDefine.hpp"
#include <cmath>
#include <string>

_NNU_DGGS_BEGIN
    namespace Base
    {
        // 常量定义
        #define M_PI  3.1415926535
        #define earthRadius 6371.00877140
        #define PRECISION 0.0000000000005
        #define tol 0.000000000000001
        #define anglePrecision 0.0000000005

        #define M_2PI       6.28318530717958647692528676655900576839433
        #define M_ONE       1.00000000000000000000000000000000000000000
        #define M_ZERO      0.00000000000000000000000000000000000000000
        #define M_HALF      0.50000000000000000000000000000000000000000
        #define M_SIN60     0.8660254037844386467637231707529361834714

        // 数学常量
        #define degToRad 0.0174532925199432957692369076848861271111
        #define radToDeg 57.29577951308232087679815481410517033240547
        #define doublePI 6.28318530717958647692528676655900576839433
        #define PIHalf (M_PI / 2)
        #define sqrtThree 1.7320508075688772935274463415058723669428
        #define sqrtThreeHalf 0.8660254037844386467637231707529361834714

        // 内切球和六边形参数
        #define offsetY (sqrtThree / 4.0)
        #define offsetX 0.25

        #define offsetYDouble (offsetY * 2.0)
        #define skewFactor (sqrtThree / 3.0)

        // 施耐德投影参数
        #define r1Value 0.9103832815
        #define r1Squared (r1Value * r1Value)
        #define dhAngle (37.37736814 * degToRad)
        #define ghAngle (36.0 * degToRad)

        #define GH            0.6283185307179586231995926937088370323181

        #define cot30        1.7320508075688774152212090484681539237499
        #define tanDh        0.7639320224822535587705374382494483143091
        #define cosDh        0.7946544722986497166772323907935060560703
        #define sinGh        0.5877852522924731371034567928290925920010
        #define cosGh        0.8090169943749474512628694355953484773636
        #define originXOffset 0.6022955029
        #define originYOffset 0.3477354707
        #define icosaEdgeLength (2.0 * originXOffset)

        // 菱形格网参数
        #define rhoyOffset offsetY
        #define rhoyOffsetDouble (rhoyOffset * 2.0)
        #define rhoxOffset offsetX
        #define rhoxOffsetComplement (1.0 - rhoxOffset)

        // 无效数据值
        #define  noData (-9999)


        /// <summary>
        /// 三维向量
        /// </summary>
        struct Vec3D
        {
            long double x;
            long double y;
            long double z;
        } ;

        /// <summary>
        /// 经纬度高程坐标
        /// </summary>
        struct GeoCoord
        {
            double lon;
            double lat;
            double elev = noData;
        };

        /// <summary>
        /// 二维点/向量（ISEA平面）
        /// </summary>
        struct Pt2D
        {
            double x;
            double y;
        };

        /// <summary>
        /// 二维点/向量（行列平面）
        /// </summary>
        struct RCPt2D
        {
            long long int i;
            long long int j;
        };


        /// </summary>
        /// 空间范围下的AABB包围盒
        /// </summary>
        struct DGGSSpaceBBox
        {
            double latMax=0;
            double latMin=0;
            double lonMin=0;
            double lonMax=0;
        };

    }

/// <summary>
/// 经纬度与高程转XYZ坐标
/// </summary>
/// <param name="lon">经度</param>
/// <param name="lat">纬度</param>
/// <param name="elev">高程</param>
/// <param name="x">x坐标</param>
/// <param name="y">y坐标</param>
/// <param name="z">z坐标</param>
    inline void LLE2XYZ(double lon, double lat, double elev, long double &x, long double &y, long double &z) {
        x = cosl(lat) * cosl(lon);
        y = cosl(lat) * sinl(lon);
        z = sinl(lat);


        if (fabsl(x) < tol) x = 0.0L;
        if (fabsl(y) < tol) y = 0.0L;
        if (fabsl(z) < tol) z = 0.0L;
    }

/// <summary>
/// XYZ坐标转经纬度与高程
/// </summary>
/// <param name="x">x坐标</param>
/// <param name="y">y坐标</param>
/// <param name="z">z坐标</param>
/// <param name="lon">经度</param>
/// <param name="lat">纬度</param>
/// <param name="elev">高程</param>
    inline void XYZ2LLE(double x, double y, double z, double &lon, double &lat, double &elev) {
        elev = 0;
        if (fabsl(z) - 1.0L < PRECISION)
        {
            if (z > 1.0L) z = 1.0L;
            if (z < -1.0L) z = -1.0L;
            lat = asinl(z);
            if ((lat == PIHalf) || (lat == -PIHalf)) lon = 0.0L;
            else lon = atan2l(y, x);
        }
        else
        {
            lon = -999;
            lat = -999;
        }

    }


/// <summary>
/// 计算两点平面距离
/// </summary>
/// <param name="p1">起点</param>
/// <param name="p2">终点</param>
/// <returns>距离</returns>
    inline long double distance(const Base::Pt2D &p1, const Base::Pt2D &p2) {
        return sqrt(pow(p2.x - p1.x, 2) + pow(p2.y - p1.y, 2));
    }

/// <summary>
/// 计算点到线的平面距离
/// </summary>
/// <param name="p">目标点</param>
/// <param name="l1">线起点</param>
/// <param name="l2">线终点</param>
/// <returns>距离</returns>
    inline long double distance(Base::Pt2D p, Base::Pt2D l1, Base::Pt2D l2)
    {
        long double A = p.x - l1.x;
        long double B = p.y - l1.y;
        long double C = l2.x - l1.x;
        long double D = l2.y - l1.y;

        long double dot = A * C + B * D;
        long double len_sq = C * C + D * D;
        long double param = dot / len_sq;

        long double xx, yy;

        if (param < 0) {
            xx = l1.x;
            yy = l1.y;
        } else if (param > 1) {
            xx = l2.x;
            yy = l2.y;
        } else {
            xx = l1.x + param * C;
            yy = l1.y + param * D;
        }

        long double dx = p.x - xx;
        long double dy = p.y - yy;

        return std::sqrt(dx * dx + dy * dy);
    }

/// <summary>
/// 计算向量点积
/// </summary>
/// <param name="vec1">向量一</param>
/// <param name="vec2">向量二</param>
/// <returns>点积</returns>
    inline long double dotProduct(const Base::Pt2D &vec1, const Base::Pt2D &vec2) {
        return vec1.x * vec2.x + vec1.y * vec2.y;
    }

/// <summary>
/// 计算向量的叉积
/// </summary>
/// <param name="vec1">向量一</param>
/// <param name="vec2">向量二</param>
/// <returns>叉积</returns>
    inline long double crossProduct(const Base::Pt2D &p1, const Base::Pt2D &p2) {
        return p1.x * p2.y - p1.y * p2.x;
    }

/// <summary>
/// 计算向量的叉积
/// </summary>
/// <param name="origin">原点</param>
/// <param name="end1">向量一终点</param>
/// <param name="end2">向量二终点</param>
/// <returns>叉积</returns>
    inline long double crossProduct(const Base::Pt2D &origin, const Base::Pt2D &end1, const Base::Pt2D &end2) {
        return (end1.x - origin.x) * (end2.y - origin.y) - (end1.y - origin.y) * (end2.x - origin.x);
    }

/// <summary>
/// 计算角度
/// </summary>
/// <param name="vecStart">线段起点</param>
/// <param name="vecEnd1">线段一终点</param>
/// <param name="vecEnd2">线段二终点</param>
/// <returns>角度</returns>
    inline double calAngle(const Base::Pt2D &endPoint1, const Base::Pt2D &startPoint, const Base::Pt2D &endPoint2) {
        //计算两线段与y轴夹角
        long double dx1 = startPoint.x - endPoint1.x;
        long double dy1 = startPoint.y - endPoint1.y;
        long double dx2 = endPoint2.x - endPoint1.x;
        long double dy2 = endPoint2.y - endPoint1.y;
        long double angle1 = atan2(dy1, dx1);
        long double angle2 = atan2(dy2, dx2);
        long double result;

        //保证计算angle1和angle2在同一侧
        if (angle1 * angle2 >= 0)
            result = abs(angle1 - angle2);
        else
            result = abs(angle1) + abs(angle2);

        //保证计算角度在0-180度
        if (result > M_PI)
            result = M_PI * 2 - result;
        return (double) result;
    }


/// <summary>
/// 判断线线是否相交
/// </summary>
/// <param name="a1">线段一起点</param>
/// <param name="a2">线段一终点</param>
/// <param name="b1">线段二起点</param>
/// <param name="b2">线段二终点</param>
/// <returns></returns>
    inline bool Intersect(Base::GeoCoord a1, Base::GeoCoord a2, Base::GeoCoord b1, Base::GeoCoord b2) {
        //快速排斥实验
        if ((a1.lon > a2.lon ? a1.lon : a2.lon) < (b1.lon < b2.lon ? b1.lon : b2.lon) ||
            (a1.lat > a2.lat ? a1.lat : a2.lat) < (b1.lat < b2.lat ? b1.lat : b2.lat) ||
            (b1.lon > b2.lon ? b1.lon : b2.lon) < (a1.lon < a2.lon ? a1.lon : a2.lon) ||
            (b1.lat > b2.lat ? b1.lat : b2.lat) < (a1.lat < a2.lat ? a1.lat : a2.lat)) {
            return false;
        }
        //跨立实验
        if ((((a1.lon - b1.lon) * (b2.lat - b1.lat) - (a1.lat - b1.lat) * (b2.lon - b1.lon)) *
             ((a2.lon - b1.lon) * (b2.lat - b1.lat) - (a2.lat - b1.lat) * (b2.lon - b1.lon))) > 0.00001 ||
            (((b1.lon - a1.lon) * (a2.lat - a1.lat) - (b1.lat - a1.lat) * (a2.lon - a1.lon)) *
             ((b2.lon - a1.lon) * (a2.lat - a1.lat) - (b2.lat - a1.lat) * (a2.lon - a1.lon))) > 0.00001) {
            return false;
        }
        return true;
    }

/// <summary>
/// 判断线面是否相交
/// </summary>
/// <param name="polygon">多边形</param>
/// <param name="pointNum">多边形点数</param>
/// <param name="point1">线段起点</param>
/// <param name="point2">线段终点</param>
/// <returns></returns>
    inline bool Intersect(Base::GeoCoord *polygon, int pointNum, Base::GeoCoord point1, Base::GeoCoord point2) {
        for (int i = 0; i < pointNum; i++) {
            Base::GeoCoord point3 = polygon[i];
            Base::GeoCoord point4 = polygon[(i + 1) % pointNum];
            if (Intersect(point3, point4, point1, point2))
                return true;
        }
        return false;
    }


/// <summary>
/// 判断面与面是否相交
/// </summary>
/// <param name="srcPolygon">原多边形</param>
/// <param name="srcPointNum">原多边形点数</param>
/// <param name="tagPolygon">目标多边形</param>
/// <param name="tagPointNum">目标多边形点数</param>
/// <returns></returns>
    inline bool Intersect(Base::GeoCoord *srcPolygon, int srcPointNum, Base::GeoCoord *tagPolygon, int tagPointNum) {
        int pointCount = tagPointNum;
        for (int i = 0; i < pointCount; i++) {
            Base::GeoCoord point1 = tagPolygon[i];
            Base::GeoCoord point2 = tagPolygon[(i + 1) % pointCount];
            if (Intersect(srcPolygon, srcPointNum, point1, point2))
                return true;
        }
        return false;
    }

/// <summary>
/// 判断点是否在面内
/// </summary>
/// <param name="polygon">多边形</param>
/// <param name="pointNum">多边形点数</param>
/// <param name="point">点</param>
/// <returns></returns>
    inline bool isPointInPolygon(Base::GeoCoord *polygon, int pointNum, Base::GeoCoord point) {
        Base::GeoCoord leftBoundPoint = {-180.0, point.lat};
        int segmentCount = 0;
        for (int i = 0; i < pointNum - 1; i++) {
            if (Intersect(polygon[i], polygon[i + 1], leftBoundPoint, point))
                segmentCount++;
        }
        if (segmentCount % 2 == 0)
            return false;
        else
            return true;
    }


_NNU_DGGS_END

